#include "Random.hpp"
#include <armadillo>

#include "Kernel_IMH.hpp"
#include "KernelMALA.hpp"
//#include "Misc/Laplace.hpp"
#ifndef _NMOVE
#define _NMOVE 1
#endif
#define _NUk 2
#define SCALE 0.01//(11/9)
#define TARGET 0.23
#define RC 2.38
using namespace arma; 

template<class Density, class Proposal>//Proposal avec r et d
class Kernel_MH{
	public:
		Kernel_MH(Density *d, Proposal *p)
		{
			_d=d;
			//_p=malloc(sizeof(Proposal));
			_p=p;
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_Sig=1;
			
		}
		~Kernel_MH(){
			delete runif;
		}
		Kernel_MH<Density,Proposal>(const Kernel_MH<Density,Proposal>& X ){
			cout << "Hello";
		}
	
/*		Kernel_MH<Density,Proposal>& Kernel_MH<Density,Proposal>::operator=(const Kernel_MH<Density,Proposal>& X){
			cout << "Hello";	
		}
*/
		void Move(mat *X,int _n)//Move matrix convert to multidimentional theta 
		{		//, moves each column dimentionis in line
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			//if(_bool==1)
		//	{
	//			Moved.zeros();
//			}
			mat Y(n,m);
			Y.fill(0);
			mat xt;
			_p->Set_n(m);
			_sum=0;
			_i=0;
			for(int i=0;i<n;i++)
			{
				Y(i,span::all)=(_p->r(i)).t()+(*X)(i,span::all);//Y[i]=p->r(X[i])
				mat xt=Y.row(i);
				xt=Accept(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;				
				_i++;
			}
			_accept=(double)_sum/n;

			cout << "accept = " << (double)_sum/n;

			//////////
			//delete[] y;
		}
		void testMove(mat *X,int _n)//Move matrix convert to multidimentional theta 
		{		//, moves each column dimentionis in line
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			if(_bool==1){
				Moved.zeros();
			}
			mat Y(n,m);
			Y.fill(0);
			mat xt;
			_p->Set_n(m);
			_sum=0;
			_i=0;
			Distribution::InvKhi2 Rkhi(1,_NUk);
			_Sig=1;
			_p->Set_s(_Sig*_s);
			mat khi(n,1);
			for(int i=0;i<n;i++)
			{
				khi(i,0)=sqrt(SCALE*_NUk*(Rkhi.r(1))(0,0));
				_p->Set_s(khi(i,0)*(RC/m)*_s);
				Y(i,span::all)=(_p->r(i)).t()+(*X)(i,span::all);//Y[i]=p->r(X[i])
				mat xt=Y.row(i);
				xt=Accept(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;				
				_i++;
			}
			//cout << khi;
			double *Yt=new double[n];
			mat2vector<double>(Moved,Yt);

			
			_accept=(double)_sum/n;
			cout << "accept = " << (double)_sum/n;
			if(_accept>0 &_accept<1)
			{
				mat mu(2,1);
				mu.zeros();

				mat bar(2,2);
				bar.fill(50);
				mat S=diagmat(bar);
				mat khi2=add_one(log(khi),1);
			//	cout << khi;	
			//	cout << "mean" << mean(khi,0); 
				boost::math::normal_distribution<> dn(0,1); 

				Laplace fL;
				try{
					fL(khi2,Yt,mu,S(0,0),mu);
					mat param=fL.Get_Mu();
					_Sig=(2.38/m)*exp((quantile(dn,TARGET)-param(1,0))/param(0,0));
				}catch(std::exception e)
				{
					cout << "Probit regression not possible";
					_Sig=2.38/sqrt(m);				
				}
			}else{
				_Sig=2.38/sqrt(m);				
			}
				_p->Set_s(_Sig*_s);
			
			//cout << khi2;
		 	//cout <<"quantile "<<quantile(dn,TARGET);
			delete Yt;

			//Set_Sig here
			//////////
			//delete[] y;
		}
		mat Accept(mat X, mat Y,int n){
	
			mat res=X;
			if(_bool==1)
			{
				double rho=0;
				double temp=_d->Density_utn(Y.t(),n);
				double temp2=(_d->Density_utn(X.t(),n));
				rho=(temp)-temp2;
				double u=log((*runif)());
				if(u<rho)
				{
					Dens(_i,0)=temp;
					res=Y;
					_sum+=1;
					Moved(_i,0)+=1;
				}else{
					Dens(_i,0)=temp2;
				}
			}else{
				
				double rho=0;
				double temp=_d->Density_utn(Y.t(),n);
				rho=temp-(Dens(_i,0));
				double u=log((*runif)());
				if(u<rho)
				{
					res=Y;
					_sum+=1;
					Dens(_i,0)=temp;
					Moved(_i,0)+=1;
				}

			}
			return res;
		}
		void Set_s(mat X){	
			mat s1=_p->Get_s();
			_s=0.5*(cov(X)+s1);
			_p->Set_s(_Sig*_s);
		}
		double Get_accept(void){return _accept;	}
	 	void Set_Sigma(double x){_Sig=_Sig*x;}
		void Set_bool(int b){_bool=b;}
		void Set_n(int b){_n=b;
			mat De(_n,1);
			mat M(_n,1);
			Moved=M;
			Dens=De;
		}
		double nMove(void)
		{
			double b=1;
			for(int i=0;i<_n;i++)
			{
				if(Moved(i,0)<_NMOVE)
				{
					b=0;
				}
			}
			return b;
		}

		double nMove2(mat X)
		{
			double b=1;
			for(int i=0;i<_n;i++)
			{
				if(Moved(i,0)<_NMOVE)
				{
					b=0;
				}
			}
			return b;
		}
	private:
		int _sum;
		Density *_d;
		Proposal *_p;
		int _bool;
		RandomG::Random<boost::random::uniform_01<> > *runif;
		mat Dens; 
		int _n;
		double _accept;
		double _Sig;
		int _i;
		mat Moved;
		mat _s;
};


